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Abstract 

We present a framework to transform the problem of finding a Lyapunov function of a Chemical 
Reaction Network (CRN) in concentration coordinates with arbitrary monotone kinetics into finding a 
common Lyapunov function for a linear parameter varying system in reaction coordinates. Alternative 
formulations of the proposed Lyapunov function is presented also. This is applied to reinterpret previous 
results by the authors on Piecewise Linear in Rates Lyapunov functions, and to establish a link with 
contraction analysis. Persistence and uniqueness of equilibria are discussed also. 
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1 Introduction 


Chemical Reaction Networks Theory (CRNs) is a multi-disciplinary area of research connecting engineering, 
mathematics, physics and systems biology. Despite diverse applications in engineering and science, the recent 
interest in CRNs is mainly due to the emergence of systems biology. 

In the context of systems biology, CRNs are key to understand complex biological systems at the cellular 
level by explicitly taking into account the sophisticated network of chemical interactions that regulate cell 
life. This is because all major biochemical networks such as signalling pathways, gene-regulatory networks, 
and metabolic networks are naturally cast in the framework of CRNs. 

A major obstacle in biochemical networks is the very large degree of uncertainty inherent in their modeling 
as well as the variability or poor knowledge of the parameters involved. Thus, a useful analysis of these 
networks shall be robust, for instance with respect to arbitrary variations in the values of parameters. In 
other words, it is desirable to provide nontrivial conclusions based solely on the structural properties of the 
network. However, it is natural to expect that this is not always possible, since dynamics mw be subject 
to bifurcations which are entirely dependent on parameter values. Fortunately, many results [l|, Q, 0] 
identified wide classes of CRNs in which such analysis is indeed possible, where it has been shown that 
some major dynamical properties such as stability, persistence, monotonicity, etc can be determined based 
on structural information only, and regardless of the parameters involved. 


I 


Early work [1|, [3j studied Lyapunov stability for weakly reversible networks which meet a certain graphical 
condition known as the zero deficiency property. It was shown that for this class of networks with mass- 
action kinetics, there exists a unique equilibrium in the interior of each invariant manifold, which is locally 
asymptotically stable regardless of the constants involved. In addition, the absence of boundary equilibria 


implies global asymptotic stability [^. Nevertheless, whether this is true without this additional assumption 
is still an open problem. 

A more recent approach for the stability problem utilized monotonicity Q . If a derived system from the CRN 
is shown to be monotone, stability theorems for mono_tone _^^ems can be used. Piecewise linear Lyapunov 
functions based methods has been proposed recently 


Hflili- 


In a previous work [7j, l8( , the authors proposed a direct approach to the problem, where Piecewise-Linear 
in Rates (PWLR) Lyapunov functions were introduced. In addition to their simple structure, these functions 
are robust with respect to kinetic constants, and require mild assumptions on the kinetic rates, mass-action 
being a special case. 

In this paper, we generalize this approach by transforming the problem of finding a Lyapunov function 
expressed in concentration coordinates for networks with arbitrary monotone kinetics into finding a common 
Lyapunov function for a linear parameter varying system written in reaction coordinates. Furthermore, sev¬ 
eral related Lyapunov functions are presented. As a result, we link the PWLR L yap unov functions introduced 


in [7( with results known in literature for piecewise linear Lyapunov functions [IflI . I 111 . Il2 1. Furthermore, we 
show to interpret our results in terms of contraction analysis and variational dynamics. Related results on 
persistence and uniqueness of equilibria are also presented. Note that some of the results in this paper were 


Ill- 


published in 

This paper is organized as follows. In Section 2, we present the background and assumptions. Section 
3 defines what we mean by a Robust Lyapunov function, and introduces the uncertain systems framework. 
PWLR Lyapunov functions are introduced in Section 4, uniqueness of equilibria and persistence is presented 
in Section 5. Section 6 discuss the relationship with contraction analysis. The proofs are presented in the 
Appendix. 


Notation Let A C M" be a set, then A°,A,3A,coA denote its interior, closure, boundary, and convex 
hull, respectively. Let x G M" be a vector, then its £oo-norm is ||a;||oo = maxi<i<„ |a;i|. The inequalities 
a:>0, a;>0, a::l:g>0 denote elementwise nonnegativity, elementwise nonnegativity with at least one positive 
element, and elementwise positivity, respectively. A is a signature matrix if A G {±1}"^" and A is diagonal. 
Let A G then ker(A) denotes the kernel or null-space of A, while Im(A) denotes the image space 

of A. A G is Metzler if all off-diagonal elements are nonnegative. The set of n x n real symmetric 

matrices is denoted by S". Let A G then A > (>)0 denotes A being positive semi-definite (definite), 
respectively. A ^ 0 denotes A being nonnegative. The all-ones vector is denoted by 1 , where its dimension 
can be inferred from the context. Let C its conic hull denotes the set KXi : Xi G R+}. 

Let V : D —)■ R, then the kernel of V is ker(F) = F“^(0). TM is the tangent bundle of the differentiable 
manifold M. 


2 Background on Complex Reaction Networks 


We summarize in this section main definitions and notions of CRNs 
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2.1 Ordinary Differential Equations Formulation 

A Complex (or Chemical) Reaction Network (CRN) is defined by a set of species y = {Xi, A„}, and set 
of reactions = {Ri,R^}. Each reaction is denoted as: 


R, 


Si=l OHjXi > PijXi, j — 1, I/, 


( 1 ) 


where , Pij are nonnegative integers called stoichiometry coefficients. The expression on the left-hand side 
is called the reactant complex, while the one on the right-hand side is called the product complex. The forward 
arrow refers to the idea that the transformation of reactant into product is only occurring in the direction 
of the arrow. If it occurs in the opposite direction also, the reaction is said to be reversible and its inverse 
is listed as a separate reaction. The reactant or product complex can be empty, though not simultaneously, 
which is used to model external inflows and outflows of the CRN. 

A nonnegative concentration Xi is associated to each species Xi. Each chemical reaction R^ is specified 
by a rate function Rj : R" —R+ which is assumed to satisfy the following: 

Al. it is a function, i.e. continuously differentiable; 


A2. Xi = 0 ^ Rji^) = Oj for all i and j such that aij > 0; 


A3, it is nondecreasing with respect to its reactants, i.e 

dRj , , J > 0 : a^j > 0 

^ I = 0 : = 0 

A4. The inequality in ([2]) holds strictly for all x G R". 


( 2 ) 


A widely-used expression for the reaction rate function is the Mass-Action which is given by: Rj{x) = 
kj nr=i > with the convention 0° = 1. 

The stoichiometry coefficients are arranged in an n x ^ matries A = [aij], B = [jdij]. They can be 
subtracted to yield T = [71 .. 7 n]^ which is called the stoichiometry matrix, which is defined element-wise as: 
[rjij = Pij — CHij. 

Therefore, the dynamics of a CRN with n species and v reactions is described by a system of ordinary 
differential equations (ODEs) as: 


x{t) = rR{x{t)), Xo := a:(0) € R" (3) 

where x{t) is the concentration vector evolving in the nonnegative orthant M" , T G R"^*^ is the stoichiometry 
matrix, R{x{t)) G Rljl is the reaction rates vector. 

Note that (l3|) belongs to the class of nonnegative systems, i.e, R" is forward invariant. In addition, the 
manifold “i^xo •= ({^0} + Im(r)) n R" is forward invariant, and it is called the stoichiometric compatibility 
class associated with Xo. 

Furthermore, the graph is assumed to satisfy: 


AG. There exists v G kerb such that u ^ 0. This condition is necessary for the existence of equilibria in 
the relative interior of stoichiometric compatibility classes. 


The set of reaction rate functions, i.e. kinetics, satisfying A1-A4 for a given A, B is denoted by J^fA- A CRN 
family yA,B is the triple {y,t%,yA). 
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2.2 Graphical Representation 

A CRN can be represented via a bipartite weighted directed graph given by the quadruple (Vs, Vr, E, W), 
where Vs is a set of nodes associated with species, and Vr is associated with reactions. 

The edge set E d V x V \s defined as follows. Whenever a certain reaction Rj is given by o, then 
(Aj, Rj) G E for all Xi’s such that > 0. That is, {Xi, Rj) G A iff aij > 0, and we say in this case that Rj 
is an output reaction for Xi. Similarly, we draw an edge from Rj G Vr to every Xi G Vs such that /3ij > 0. 
That is, (Rj , Ai) G E whenever Pij > 0, and we say in this case that Rj is an input reaction for Xi. Note 
that there are no edges connecting two reactions or two species. The weight function W : i? —>■ N assigns to 
each edge a positive integer as W{Xi,Kj) = aij, and W{Rj,Xi) = Pij. Hence, the stoichiometry matrix T 
becomes the incidence matrix of the graph. 

The set of output reactions of P is denoted by A(P). A nonempty set P C Vs is called a siphon [? ] 
if each input reaction associated to P is also an output reaction associated to P. A siphon is a deadlock if 
A(P) = Vr. a siphon or a deadlock is said to be critical if it does not contain a set of species corresponding 
to the support of a conservation law. 

3 Robust Lyapunov Functions and Linear Inclusions 

3.1 Robust Lyapunov Functions 

In order for the stability analysis of CRNs to be independent of the specific kinetics, we aim at constructing 
Lyapunov functions which are dependent only on the graphical structure, and hence are valid for all reaction 
rate functions that belongs to Jttr- Therefore, we state the following definition: 

Definition 1 (Robust Lyapunov Function). Given ([3]). Let V : ^ 1R+ be locally Lipschitz, and let 

—>■ be a function, where Xg is an equilibrium. Then, (V, Wn^xe) said to induce a Robust 

Lyapunov Function (RLF) with respect to the network family PyA.B if for any choice of R G G M", 

the function Vh,x„ = V o Wr^x^ is, 

1. Positive-Definite: V{x) > 0, and V{x) = 0 if and only if R{x) G kerT. 

2. Nonincreasing: V(x) < 0 for all x G ’^xe ■ 

A network for which an RLF exists is terms a Graphically Stable Network (GSN). 

Remark 1. As will be seen afterwards, the function V used in the definition of the Lyapunov function 
(through composition) is invariant with respect to the specific network realization in JFa, while the function 
Wr^xs is allowed to depend on the kinetics of the network. Two main examples of the function Wr^x^ 
Wr^x^x) = R(x), and Wr^x,,{x) = x — Xe, and they will be used afterwards. 

As an abuse of notation, we will call the parametrized Lyapunov function Vr^x^ RLF. 

Remark 2. The time-derivative in the definition above is the upper right Dini’s derivative H/; 

(4) 

which is finite for all x since V is locally Lipschitz. 

Since the RLF defined above is not strict, we need the following definition: 
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Definition 2. A n RLF for ^a,b is said to satisfy the LaSalle’s condition if for any choice R € J^a 
the following statement holds: 

If a solution ip(t;Xo) of ([3]) satisfies (p{t;xo) € kerV fl ^ > 0, then this implies that tp(t;xo) G Ex„ for 
all t > 0, where Ex^ C “^xo the set of equilibria for 


The following theorem adapts Lyapunov’s second method 



to our context. 


Theorem 1 (Lyapunov’s Second Method). Given ([3]) with initial condition Xo G R”, and let '^xo os the 
associated stoichiometric compatibility class. Assume there exists an RLF Lyapunov function, and suppose 
that x{t) is bounded, 


1. then the equilibrium set Ex,^ is Lyapunov stable. 

2. If in addition, V satisfies the LaSalle’s Condition, then x{f) -A- Ex„ ast^oo (i.e., the point to set dis¬ 
tance of x(t) to Ex,, tends to 0). Furthermore, any isolated equilibrium relative to '^xo is asymptotically 
stable. 

3. If V satisfies the LaSalle’s condition, and all the trajectories are bounded, then: If there exists x* € 
Ex,, which is isolated relative to '^xo then it is unique, i.e., Ex, = {x*}. Furthermore, it is globally 
asymptotically stable equilibrium relative to '^x, ■ 


Remark 3. Note that the RLF Lyapunov function considered can not be used to establish boundedness, as 
it may fail to be proper. Therefore, we need to resort to other methods so that the boundedness of solutions 
can be guaranteed. For instance, if the network is conservative, i.e the exists w G M" such that w'^T = 0, 
which ensures the compactness of ‘^x, ■ 


3.2 Uncertain Systems Framework: Reaction Coordinates 

In this subsection, the function Wr^x, is assigned to be R, as in the case of PWLR functions. As arbitrary 
monotone kinetics are allowed in our formulation of the CRN family aVa^b-, the system ([3]) with kinetics JFa 
can be viewed as an uncertain system. However, this system does not fit directly to the traditional types of 
uncertainties (e.g parameter uncertainty) known in the literature. In this subsection, we show that shifting 
the analysis of the system to reaction coordinates enables us to view it as a linear parameter varying (LPV) 
system where the existence of a common Lyapunov function for the LPV system implies the existence of a 
robust Lyapunov function for the CRN. 

Let r{t) := R{x{t)), then we have: 

r{t) = ^i^mrit)=p{t)rr{t), (5) 

where p{t) := ^{xft)). We can write p{t) as a conic combination of individual partial derivatives as follows: 

no 

— {x{t))=p{t)= pMt)Ej„ (6) 

>0 

where [p(t)]ji = Pji(t), and [Eji]j>if = 1 if {j',i') = {j,i) and zero otherwise. 

Let s denote the number of elements in the support of dR/dx, and let k : {1,.., s} -A {{i,j) : cnj > 0} be an 
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indexing map. Then, we can write (O as: 

S 

r= Pji{t)Ej,rr = (7) 

i,j-.aij>0 

where , pi{t) = Pji{t), with {i,j) = k{£), and {ej}^^^ denotes the canonical basis of Hence, 

eq. © represents a linear parameter-varying system which has s nonnegative time-varying parameters 
{pi{t),Ps{t)} where the system matrix belongs to the conic hull of the set of rank-one matrices {T^,T®}. 
Hence, we have the following definition: 

Definition 3 (Common Lyapunov Function). A function V : Rlji —> K+ is said to be a Lyapunov function for 
the linear system r = F^r if it is locally Lipschitz, nonnegative, has a negative semi-definite time-derivative, 
and ker V C ker F^ . Furthermore, V is said to he a common Lyapunov function for the set of linear 
systems {r = T^r, ...,r = r®r} if it is a Lyapunov function for each of them, and ker V = riLikerrC 

Hence, we are ready to state the main result whose proof is present in the appendix: 

Theorem 2 (Common Lyapunov Function). Given the system ([^. There exists a common Lyapunov 
function V : R+ for the set of linear systems {f = F^r,r = F®r}, if and only if (V, R) induces the 

Robust Lyapunov function parameterized as VrIx) = V{R{x)) for the CRN family ^a,b- 

Remark 4. Since the zero matrix belongs to the conic hull of {F^,...,r®}, asymptotic stability can’t be 
established by the mere existence of the common Lyapunov function. A LaSalle’s argument is needed as will 
be mentioned in the following section. 


3.3 Dual Robust Lyapunov Function: Species Coordinates 


The RLF introduced in the previous subsection was a function of R{x). We investigate now RTFs that are 
functions oi x — Xe- This will be carried out in a manner that is dual to what has been done in the previous 
subsection. 

In order to present the dual framework, an alternative representation of the system dynamics can be 
written. Given ®, and let Xe be an equilibrium. Then, there exists x"(x) such that m can written 
equivalently as: 

D D 

X = T — {x"){x - Xe),x{0) (8) 

ox 

The existence of x" := x^ F ex{x — Xe) for some Ex & [0,1] follows by applying the Mean-Value Theorem 
to R{x) around Xe- 

Let z = X — Xe, then similar to the previous section, the conic combination ([6]) can be used to rewrite ® 


as: 


OR 


z = r—{x"){x - Xe) = = y^p^(t)F. 


e=i 


=1 


(9) 


where p£{t) = -q^{x''{ x{t)), and F^ is the column of F. Therefore, the systems dynamics has been 
embedded in the linear differential inclusion with vertices ...,ri^ej^}. 

Let be a matrix whose columns are basis vectors of kerF^, the following theorem can be stated: 


Theorem 3 (Common Lyapunov Function). Given the system ([3]). There exists a common Lyapunov 
function V : R" —>■ M+ for the set of linear systems {z = (Fi^eJ^)z,..., i = (ri^e^)^}, with the constraint 
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D"^z{0) = 0 for each of them if and only if {V, Wx^), Wx^ = x — Xe induces the Robust Lyapunov function 
parameterized as Vx^(x) = V{x — Xe) for the CRN family NCa,b- 

Proof. Let V{x) = V{x — Xe), then when dVfdz exists we can write: 

■ dV dV^ , , j. ^ . fdV r 

^ ^ ^ ^ ^ pe{t) 

e=i e=i \ 

Since we assumed that F is a common Lyapunov function for the set of linear systems {z = (LijcJ^z, z = 
(Ti^ejjz} the proof can proceed in both directions in a similar way to the proof of Theorem[21 Note that 
the constraint Dz{0) = 0 is needed since Dx{t) = 0 is a constraint imposed by the structure of the original 
system ([3]) □ 

3.4 Relationship Between the Two Frameworks 

We show now that if V can be written as a certain composition, then the Lyapunov function of the form 
V(x — Xe) can be used, where Xe is an equilibrium point for ([3]). 

Hence, the following theorem can be stated: 

Theorem 4 (Relationship). Let Vi{x) = V{R{x)) be representing an RLF for the network family jPa ,B- If 
there exists V : K„ —>■ ]R_|_ such that for all r: 



V{r) = V{Tr), 


( 10 ) 


then V 2 {x) = V{x — Xe) represents an RLF for the network family pyA,B, where Xe is an equilibrium point 
for (IH). 

Proof. Condition 1 in Definition[T]is clearly satisfied. It remains to show the second condition. Let z = x—Xg. 
Then, whenever V is differentiable: 


V2(x) 


dV(x-Xe) . dV(x-Xe)^^, , 

-X = -- 1 R[x), 


dz 


dz 


Before proceeding, we state two statements that we need: First, from (fTUl) . we get [dV (r)/ dr) = {dV{Tr)/dz)T. 
Second, note that x — Xg € Im(r), hence 3R{x') such that rR{x') = x — Xg, where R{x') can always be 
chosen nonnegative by A7. Hence, when V is differentiable, we can use (|H]) to write: 


tW ^ dV{x-Xg) dR{x'') 

V2{x) = ---r—-- [X - Xg) 


dr dx 

dV(R(x')) dR{x'') 
dr dx 

dV{Rix')). 


rR{x') 


e=i 


dr 


rR{x') < 0 , 


where the last inequality follows from (l26l) . Lemma [Tbl implies that V 2 {x) < 0 for all x. □ 

Remark 5. Note that V 2 has a simpler structure that Vi since it depends on x — Xg. However, it can be 
noticed in the proof that for a specific choice of Xg, the Lyapunov function V(x — Xg) is nonincreasing only 
along trajectories that starts in '^x,. ■ 
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Remark 6. Note that if we consider the system for z = x where 


z = T^^z,Dz{Q) = 0. ( 11 ) 

then V(z) is a Lyapunov function for the above system. 

Relationship to the Extent of Reaction Formnlation 

Recall that the extent of reaction i is defined as: f{t) = f* R(x(T))dT + ^(0). If x(t) € then 3^* > 0 
such that Xe — Xo = We set ^(0) '■=£,*■ Hence, we can write: 

T^{t) = x{t) - Xe, (12) 

and 

^ = i?(xe + rC),^(0):=r, (13) 

which is the extent-of-reaction ODE representation of the dynamics of the CRN. 

Therefore, we state the following result: 

Corollary 5. Given F. Let V : Klji —>■ K+. Assume there exists V : K” —>■ M+ such that for all r, 
V{r) = H(Fr). IfV is a common Lyapunov function for {r = ...,r = Cj^jf^r} where {ii,ji) = «(€). 

Then, E('C) is nonnegative, and nonincreasing along the trajectories of ^ = R{xe + F^). 

Proof. The first two statements follow from Theorems I2I4I We prove now the third statement. Using m, 
m we get U(C) = U(a; — Xe)- Therefore, the required statement follows from the result that V{x — Xe) is 
nonincreasing along the trajectories of (jS]). 

□ 


4 Application to PWLR Lyapunov Functions 

4.1 Relationship to Previous Results 

In the previous papers 0, , the concept of Piecewise Linear in Rate (PWLR) Lyapunov functions has been 

introduced based on a direct analysis of the CRN. Such functions satisfy the conditions of Definition [U and 
hence they are Robust Lyapunov functions. In this section, we show that those results can be interpreted 
in the uncertain systems framework introduced above. This also allows us to provide alternative algorithms 
for the existence and construction of PWLR functions. 

Consider a CRN (|3]) with a F £ Two representation were discussed of the PWLR Lyapunov 

function. Given a partitioning matrix H G such that keriJ = kerP. PWLR Lyapunov functions are 

piecewise linear in rates, i.e. they have the form: V{x) = V{R{x)), where V : —>■ M is a continuous PWL 

function given as 

V{r) = \clr\, r G ±>Vfc,fc = 1 ,..,to/2, (14) 

where the regions Wk = {r G R^ : T,kHr > 0},fc = 1 ,..,to form a proper conic partition of R^, while 
are signature matrices with the property = —Urn+i-k, k = 1,..,to/ 2. The coefficient vectors of 
each linear component can be collected in a matrix C = [ci,.., cm]'^ G If the function V is convex. 



then we have the following simplified representation of V: 


V(x) = \\CR{x) 


This representation reminds of the £oo Lyapunov functions that were used for linear systems as in ll|. In 


fact, the next theorem establishes the link between the results introduced in [7|, l8j for checking candidate 
PWLR functions based on direct analysis and previous work on ^aa Lyapunov functions using the framework 
introduced in the previous section; 


Proposition 6. Given L and H. Let V = V o R be the candidate continuous nonnegative PWLR with 
C = [ci ... g Then {V,R) induces an RLF if and and only if: 


1. kerC = kerL, and 

2. there exists {R^Ye.=i C R^^t such that 


A^H = -CY, (15) 

and > 0, where . 

If V is convex, then the second condition can be replaced with, 

2) there exists {A^Yi C R™^™ Metzler matrices such that 

A^C = CY, (16) 

and A^l = 0 for all i = 1 ,.., s, where C = [C^ —C'^Y- 

Proof. The proof can be carried out by performing elementary algebraic manipulations on the results pre¬ 
sented as Theorems 4,5 in Q. The details are omitted for brevity. □ 

Remark 7. Note that by the symmetries in (di), it can he written equivalently: 

CY = ¥C, (17) 

where M is an ^ x ^ matrix which be defined by subtracting the upper ^ x ^ blocks of A^. A^ will satisfy: 


max 

k 


Xkk + |Afej I j < 0. 


(18) 




This is exactly the condition for which loo-norm Lyapunov functions need to satisfy for a linear system 
'1(1 . This shows that Theorem provides us with the framework to utilize the existing linear stability 

analysis techniques in the literature to construct robust Lyapunov functions for nonlinear systems such as 
CRNs. For example, we can verify £i Lyapunov functions of the form V{x) = ||Ci?(a;)||i directly by replacing 
condition dl) by 


max I Xkk + EiviN». 


(19) 


instead of converting them to the lao-norm form. 
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Before proceeding to the construction algorithm, we need to introduce the concept of a neighbor to 
a region. Fix k G {1,to/ 2}. Consider H: for any pair of linearly dependent rows eliminate 

hf^. Denote the resulting matrix hy H G and let Si,..,Sm the corresponding signature matrices. 

Therefore, the region can be represented as Wfe = {r|Sfe77r > 0}. The distance dr between two regions 
Wfe, Wj is dehned to be the Hamming distance between . Hence, the set of neighbors of a region Wk 

are defined as: 


A4 = {j\dr{'Wj,yVk) = 1,/ = 1, 


Equivalently, note that a neighboring region to Wk is one which differs only by the switching of one inequality. 
Denote the index of the switched inequality by the map Sfe(.) : A4 —{1, For simplicity, we use the 

notation Ske ■= Sk{i)- 

We use Theorem to show that the problem of constructing a PWLR Lyapunv function over a given 
partition, i.e. a given H, can be solved via linear programming. However, instead of encoding the nonde- 
creasingness condition into precomputed sign patterns as in Q , we use here alternative condition which are 
stated in the following proposition: 


Proposition 7. Given the system © and a partitioning matrix H G Consider the following linear 

program: 

Find Cfc,a,Cfc GM, 

k= GAfkJ= 

subject to FfkH, 

CT^ = -A^H,Xij:k > 0 , 

Ck — Cj = rjkjCTkskj hskj ^ 

6>0,l^Cfc>0,A^>0. 


Then there exists a PWLR RLF with partitioning matrix H if and only if there exist feasible solution to 
the above linear program with kerC = kerF satisfied. Furthermore, the PWLR RLF can he made convex by 
adding the constraints rjkj ^ 0. 

Remark 8. A natural choice for H is H = T. Refer to Remark 17 in Q/. 


Remark 9. The LaSalle’s condition can be verified via a graphical algorithm described in %IITF in 


4.2 The Dual PWL Lyapunov Function 

In §HI-C, it was shown that if there exists V such that V(r) = F(Fr), then there exists a dual RLF for the 
same network family. In the case of PWLR Lyapunov functions, condition 1 in Proposition [6] implies that 
this is always possible. Hence, consider a PWLR Lyapunov function defined with a partitioning matrix H 
as in d- By Proposition El and the assumption that keriL = kerP, there exists such 

that H = GF, G = BY. Similar to we can define the following regions: 

Vfe = {zl'FkGz > 0},fc = 1, ..,TO, 

where it can be seen that Vk has nonempty interior iff Wk has nonnempty interior. 

Therefore, as the pair (G, H) specify the PWLR function fully, also the pair {B, G) specifies the following 
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function: 


V{z) = b^z, whenSfeG^; > 0, 


where B = [6i, If V is convex, then it can be written in the form: Vi{x) = ||Gi?(a::)||oo. Similarly, 

the convexity of V implies that 

V 2 {x) = \\B{x - Xe)||oo, 


where the later is the Lyapunov function used in 16, 

Theorem |4] established that if V{R{x)) is an RLF, then V(x — Xe) is an RLF also. The following theorem 
shows that converse holds also for PWLR RTFs, however, it is worth recalling from Remark [S] that V{R{x)) 
is nonincreasing for all initials conditions, while V(x — Xe) is nonincreasing only on 


Theorem 8. Given ([3]). Then, if there exists G € W’^'^,B G K™ such that: 


1. (iJF, GF) defines a PWLR RLF, then {B, G) defines a dual PWL RLF. 

2. (B, G) defines a dual PWL RLF, then (BF, GF) defines a PWLR RLF. 


The proof is presented in the appendix. 

Remark 10. Note that since D"’"(x — Xe) = 0 for x G ‘^xe, then if \\B{x — Xe)||oo is a Lyapunov function, 
then ||(B -GYD'^){x - Xe)||oo is also for an arbitrary matrix Y. Furthermore, since Theorem\^ showed that 
the reaction-based and the species-based representations are equivalent; it is easier to check and construct 
RLFs in the reaction-based formulation. 


5 Uniqueness of Equilibria and Persistence 

5.1 Robust Non-singularity 

It has been shown in that the Jacobian of any network admitting a PWLR RLF is Pq, which implies that 
all principal minors are nonnegative. This can be used to show the following result: 

Theorem 9 (Robust Nonsingularity). Given (jS]). Assume that it is a graphically .stable network. If for some 
realization R G JFa, there exists an isolated equilibrium Xg relative to ^x„, the reduced Jacobian with any 
realization of the kinetics is non-singular in the interior of R” . This implies that any positive equilibrium of 
this network is isolated relative to its class. 


Proof. The proof is based on the result that the negative Jacobian is Pq for any choice oi R G JFa- Using 
the Cauchy-Binet formula 17|, let I C {1, ..,n} be an arbitrary subset so that |/| = k. The corresponding 
principle minor can be written as: 




i: 


J<Z{l,..,v},\J\=k 


det(-r/j)det(^||^^^^=^^a, 


n 


pi. 


^eLtC{i,..,s} 


where (*) refers to the fact that the sum can be expressed as a linear combination of products of pi ,..., ps. 
We claim that the coefficients are all nonnegative. To show this, assume for the sake of contradiction that 
there is some Ot, is negative. If we set all p’s to zero except the ones appearing in i**' term, then this implies 
that the corresponding principle minor can be negative; a contradiction. 
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Now, the theorem can be proven by noting that the reduced Jacobian is non-singular iff the sum of all 
k X k principle minors of the negative Jacobian is positive, where k = rank(r). Since it was assumed that 
there exists an isolated positive equilibrium, this implies that the sum of principle minors is positive for some 
choice of pi,..,pa. Since all of the principle minors are nonnegative, then at least one of them is positive. 
By AK4, that principle minor will stay positive for any choice of positive pi, ..,ps, i-e. it stays positive over 
the interior of R" . □ 


5.2 Uniqueness of Equilibria 

A main result the links injectivity of a map to the notion of P-matrices which is first shown in [l^ which 
states the a map is injective if its Jacobian is a P matrix. This implies that if an equilibrium exists, then it 
is unique. This notion has been studied extensively for reaction networks, where graphical conditions were 


provided by Banaji et al. 171, Banaji and Craciun m- 




It has been established in 


19|, Appendix B] the network can not admit multiple nondegenerate positive 


equilibria in a single stoichiometric computability class if the Jacobian is Pq. Hence, the following theorem 
follows: 


Theorem 10 (Uniqueness of Positive Equilibria of GSNs). // Ua,_b is GS, then it can not admit multiple 
nondegenerate positive equilibria in a single stoichiometric computability class. 

Furthermore, if there exists an isolated positive equilibrium Xe then it is unique relative to ■ 

Proof. Since the negative of the Jacobian is Pq, the first statement follow directly from the statement before 
the theorem. 

For the second statement. Theorem [9] implies that the existence of an isolated positive equilibrium Xe 
ensures that the reduced Jacobian is non-singular on the interior of the orthant. In order to show uniqueness, 
assume for the sake of contradiction that there exists y ^ Xe,y G such that TR{y) = 0. Then the 
fundamental theorem of calculus implies, 

/•I 3R 

Q = TR{xe)-TR{y) = T J -^{txe + {1 - t)y) {xe - y)dt = T—{x*){xe - y), 

where x* = t*Xe + {l — t*)y, and t* G (0,1). The existence of t* is implied by the integral mean-value theorem. 
Since x* G then the reduced Jacobian at x* is non-singular relative to ImT. Since Xg — y G ImT, then 
y = Xe', a contradiction. □ 


Remark 11. Note that since the Jacobian is Pq, then if arbitrary inflows and outflow were added to every 
species of a GSN, then the resulting Jacobian will be a P matrix [2^]; this is also known as the continuous- 
flow stirred tank reactor (GFSTR) version of the network. The GFSTR network will he injective. This shows 
that our framework has a direct relationship to the recent results on injectvity and P matrices of reaction 
networks \l^, \l^, [3/ . 


5.3 Persistence 


Persistence is a dynamical property that is satisfied if all trajectories that start in the interior of the positive 
orthant does not approach to the boundary. Angeli et al. 23| showed that if a conservative network does 
not have critical siphons, then it is persistent. Note that this is a graphical property, and it is independent 
of the specific realization of the kinetics involved. 
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In this section, it is shown that this can be established for conservative GSNs under some conditions 
detailed in the following theorem whose proof is presented in the appendix: 

Theorem 11 (Absence of Types of Critical Siphons). Given ([3]). Consider the network family ^a,b- 
Assume there exists a critical siphon P, and let A(P) the set of output reactions of P. Then, is a 

not GS, i.e., the network does not admit a PWLR Lyapunov function if any of the following conditions is 
satisfied: 


1. If P is a critical deadlock. 

2. If the network is conservative, and if for some realization of the network family there exists an isolated 
positive equilibrium in the interior of a proper stoichiometric compatibility class. 


3. If the network is conservative, and if kerb is one-dimensional. 


The third item and the fourth items of the previous theorem can be combined with the main result in 23j , 
and by noting the inclusion of a reverse of a reaction does not create a critical siphon we get the following 
theorem: 


Theorem 12 (Persistence of A Class of CSNs). Given ([3]). Assume the network is conservative. Then the 
network family pIa.b is persistent if 

1. kerb is one-dimensional, or 

2. if by removing the reverse of some reactions the reduced network has a one-dimensional kernal and is 
GS, or 

3. there exists an isolated positive equilibrium in the interior of some proper stoichiometric class for a 
realization of jI'a.b- 


6 Relationship to Contraction Analysis 


Contraction analysis is an approach which study stability of a trajectory with respect to other trajectories, 
and not with respect to a given equilibrium point. This area of research is old as the concept of contraction 
maps is, however, it has sparked growin g in terest in the control systems community in relationship to the 
analysis of dynamical systems 241, [^, [2^. 

One the formulation of interest to us is one which utilizes matrix measures, or logarithmic norms. 


Definition 4 (Logarithmic Norms). For a given norm ||.||* on ffi.", the associated matrix measure can be 
define as a following for a matrix A € 


fj.^,(A) := limsup 
h—>-0+ 


||/ + — 1 
h 


( 20 ) 


Note that the same definition applies if was ||.||* a semi-norm. 


Remark 12. The logarithmic norm can be evaluated for the standard norms, in fact, the following expression 
can be used for the £oo norm: 


Moo (A) = max 

i 


+ X! 


( 21 ) 
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Note that this is an identical expression to one appearing in 


For a dynamical system, negativity of the logarithmic norm can be linked to contraction. This result 
were stated in different forms and have long history, refer to the tutorial for a more details. We state 
the result as the following: 

Theorem 13 (Aminzare and Sontag [^1. Consider a dynamical system x = f{x) defined on a convex 
subset o/K". Let ||.||* he a norm on K", and assume that 

Vx./r* ( ^ I < c. 

Then, for any two solutions ip{t] x), ip(t;y) of the dynamical system, the following condition holds: 

Wififx) - (/3(t;y|U < e'=*||(/5(0;a:) - (p{0;y)\\^. (22) 


Note that if c < 0, then the solutions of the system will be exponentially contracting. If c = 0, then 
the system will be non-expansive. The choice of the norm plays a crucial rule even with respect to diagonal 
weighings. Russo et al. [27j has applied the result above to CRNs by weighting the ^i-norm with a diagonal 
matrix. 

This motivates the following question: since our convex PWLR Lyapunov functions are semi-norms 
weighted by a non-square matrix, is it possible to find a link with contraction analysis?. The link between 


the logarithmic norms and norm-based Lyapunov functions has been established before |15| . Therefore, it 
might expected that such a relationship exists with contraction analysis. Indeed, the answer is yes. The 
following theorem states precisely the relationship between our convex PWLR functions and contraction 
analysis, while the proof is presented in the appendix: 


Theorem 14 (Relationship to Contraction Analysis). 1. Given the extent of reaction representation of 

CRNs ^ = R{xe + r^)- Assume that there exists convex PWLR function F(^) = HC^Hoo; and let fic 
be the logarithmic norm associated. Let the associated Jacobian be: Ji(^) := Then, 


Ve,A^c(Ji(O)<0. 


Hence, the system is non-expansive along directions orthogonal to kerP. 

2. Given the ODE x = TR{x). Assume that there exists convex PWLR function F(^) = \\B(x — Xe)\\oo, 
and let (Xb be the logarithmie norm associated. Let the associated Jacobian he: J 2 {x) := Then, 

\/x e lJ-BiJ 2 {x)) < 0. 

Hence, the system is non-expansive in each stoichiometric class ^xe- 


6.1 Variational Dynamics and LaSalle Argnment for Contraction Analysis 

In a recent work, Forni and Sepulchre proposed a Lyapunov framework for contraction analysis using a 
so called Finsler structure. In order to minimize the background needed, we apply it directly to our context. 
The Finsler-Lyapunov function for the system m is Vp : TRf —>■ R+. If contraction analysis was carried 
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out with respect to the norm: ||.||* : ^ i—^ HC'^Hoo, then the correspoinding Finsler-Lyapunov function is: 


Vf{SO = IICJ^IIoo. 


The Finsler structure is given by the mapping 
function for the variational system: 


||C'<5Clloo- Therefore, Vp can considered as a Lyapunov 


ox 

The Finsler structure induces a distance function, which is in this case dF{x,y) = 11(7(0: — y)||oo- Hence, if Vp 
was strictly decreasing then this implies that this system is incrementally asymptotically stable with respect 
to the distance function 2^, which is equivalent to the result given by Theorem 1131 However, if the strict 


decreasingness is not available, the Finsler-Lyapunov framework for contraction analysis has the advantage 
of accommodating a LaSalle’s invariance principle that can be used to show strict contraction. In fact, the 
same algorithm given in [^, §HI-F] can be used for the Finsler-Lyapunov function. The following theorem 
states the result: 

Theorem 15 (Strict Contraction). Given ([3]). Assume that V{x) = ||(7i?(x)||oo is a PWLR Lyapunov 
function that satisfies condition i. Then, the trajectories of (1131) are exponentially contractive with respect 
to the norm ||.||* : f i— ||(7^||oo in directions orthogonal to kerF. 


Proof. As per 2^, Theorem 2], we need to show that if a trajectory lives in kerVp, then it is an equilibrium 
for the variational system above. Since Vf is a piecewise function it can be studied per partition regions as 
before. Hence let 

VFm = elk = el^^St 


Note that the expression above is analogous to [8|, eq. (25)]. In fact, the arguments of the proof of Proposition 
7 in Q can be replicated to show that C5 i implies that when S£,{t) G ker Vp for all t >0, then TSf,{t) =0. □ 


Appendix: Proofs 


Before proving the results of the paper, we need to state and prove the following Lemma: 
Lemma 16. Let x := f{x), and let V : R" R+ be a locally Lipschitz function such that: 


dV(x) ^ , dV(x) 

—-- fix) < 0, whenever —-- 

dx ^ - dx 


exists, 


then V{x) < 0 for all x. 

Proof. Since V is assumed to be locally Lipschitz, Rademacher’s Theorem implies that it is differentiable 
(i.e., gradient exists) almost everywhere [^ . Recall that for a locally Lipschitz function the Clarke’s gradient 
at X can be written as dcV{x) := codV{x), where: 


dV{x) := {p e R" : 3xi —>■ xwii\\dV{xi)/dx exists, such that 

p = lim dV[xi)/dx\ . 
i—¥OQ J 

Let p G dV(x) and let be the corresponding sequence. By the assumption stated in the lemma, 

{dV{xi)/dx)f{xi) < 0, for all i. Hence, the definition of p implies that p'^f{x) < 0. Since p was arbitrary. 
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the inequality holds for all p £ dV{x). 

Now, let p £ dV(x) where p = XiPi is a convex combination of any pi, ...,pn+i £ dV(x). By the inequality 
above, f(x) = Xi{pf f{x)) < 0. Hence, p'^f{x) < 0 for all p £ dV{x). 

As in [29j, the Clarke’s derivative of H at cc in the direction of f{x) can be written as D'^f^^^Vix) = 
max{p^/(a;) : p £ dV{x)}. By the above inequality, we get ^ 0 for all x. Since the Dini’s 

derivative is upper bounded by the Clarke’s derivative, we finally get: 


V{x) := limsup 

/t^0+ 


V{x + hf{x)) — V (x) 


< limsup 

J-0 + 
y->-x 


V(y + hf(x)) -V{y) 


=: D 




for all X. 


□ 


Proof of Theorem We show the existence of the common Lyapunov function implies the existence of 
the RLF. Nonnegativity of V follows from the nonnegativity of V. Let {i,j) = k{^), recall that L^ = ejjf, 
hence kerC = = kerb. Therefore, R{x) £ kerb iff ri?(a:) = 0, which establishes the positive- 

definiteness of V. 

We assumed that V has a negative semi-definite time-derivative for every linear system in the considered 
set. Hence, when V is differentiable, we can write {dV/dr)T^r < 0, £ = 1, s. Hence, for any p ^,£ R+: 


^ QY 

p^-—T^r < 0, when(9H(r)/9r) exists. 
„ ur 
e=i 


Therefore, when V is differentiable: 


’'<*> = as & 


OR, 



ip-.OLij >0 


dxi 


{x)Eji TR{x), 


(23) 


(24) 


where dV/dR := 
Now, denote p^ = 


(dV/dr 


^^(x), which is nonnegative by A3. This allows us to write: 


r—R{x) 


v{x) = Y,p^—e,j:r{x) 

1^1 

^ ^ QY 

= < 0, for almost alia:. 


(25) 

(26) 


The last inequality follows from (l23l) . Using Lemma fTHl V{x) <0 for all x, and for all R £ JEt- 

In order to show the other direction, almost all properties outlined in Definition [3] are clearly satisfied, 
we just show nonincreasingness. Assume that there exists £ such that V (r) is not nonincreasing along the 
trajectories of r = L^r. Consider the corresponding term in pK]) . Since V{R{x)) is a Lyapunov function for 
any choice of admissible rate reaction function R, choose p^ = to be large enough such that V{x) > 0 
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for some x; a contradiction. 


Proof of Theorem [8] The first statement follows from Theorem^ In order to show the second statement, 
let V 2 (x) = b^(x — Xe), for x — Xe € Vk- We will show that Vi{x) = cJ^R{x), for R{x) G Wfe is nondecreasing. 
Without loss of generality, the partition matrix can be written in the form: G = [I G^]^. This representation 
implies that the sign of a: — a;e is determined in every region Vk- Now, assume that x — Xe & V^, then: 

^ 2 ( 3 ^) =b^TR{x) = clR{x)<0 = clR{xe), for alii? G J^a- 

Let Rj{x) G suppi?(a:), and let > 0. Since R is nondecreasing by A3, if sgn(a;i — XeJ sgn(cfcj) > 0, there 
exists R G JCi such that V 2 {x) > 0. Hence, this implies that the inequality sgn(cfcj) sgn(a;i — a^e*) < 0 holds. 
Fix j, if there exists zi ,*2 such that ai^jjUi^j > 0 and sgn(a;ij — ) sgn(a;i 2 — ) < 0, then akj ■= 0. 

Otherwise, akj ■= sgn(xi — XgJ for some i such that > 0. 

Hence, in order to have V 2 {x) < 0 for all R G JVa we need that akj{xi — XgJ >0 whenever x — Xg G Vfc, for 
all fc, j, i with aij > 0. By Farkas Lemma [^, this is equivelant to the existence of Xkji G K+, Qkji G RL 

akjef = Xl^,EkG + Clj,D, (27) 

where G is a matrix whose columns are basis vectors for kerF^. 

If we multiply both sides of (E7)) by F from the left, then we get condition C4 in [^, Theorem 4] which 
necessary and sufficient for Hi(x) = ^(c^i?(x)) <0. ■ 

Proof of Theorem [11] Assume P is a critical siphon for the petri-net associated with F, and let Up = |P|. 
Let A(P) be the set of output reactions of P, and let z/p = |A(P)|. 

Before we prove item 1 of Theorem II11 the following lemma is needed. 

Lemma 17. Consider a network family ^a,b- LC P be a set of species that does not contain the support of 
a conservation law; let its indices be numbered as {!,...,rip}. Then, there exists a nonempty-interior region 
{rlEjjFr > 0} with a signature matrix T,k that satisfies aki = ■■■ = aknp = 1- 

Proof. Assume the contrary. This implies that ri^Ti{R\jJ'R > 0} Q n(L„p^i{i?|tTi 7 j^i? > 0} = 0 for all 
possible choices of signs ai = ±1. However, R’’ can be partitioned into a union of all possible half-spaces of 
the form R > 0}. Therefore, this implies that R > 0} = 0. By Farkas Lemma, 

this implies that there exists A G R* satisfying A > 0 such that [A^OjF = 0. Therefore, P contains the 
support of the conservation law [A^ 0]^; a contradiction. □ 

Therefore, we can state the proof of the first item: 

Proof of Theorem MlV l). Without loss of generality, let {!,...,np} be the indices of the species in P. 
Using Lemma 1171 there exists a nonempty-interior sign region Sk, 1 < fc < mg with a signature matrix 
that satisfies Oki = ■■■ = aknp = 1- Since A(P) = this implies bkj < 0 for all j = 1, ..,n. However, this 
is not allowable by [^, Theorem 9] since (f[BkV < 0 for all v G kerb n R" and for any choice of admissible 

Cfe. □ 

In order to proceed, an existence result of equilibria is needed: 
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Lemma 18. Consider a network family ^a,b ■ Let P be a critical siphon, and let 'Pp fee the associated face. 
If the network is conservative, then for any proper stoichiometric compatibility , there exists an equilibrium 
Xe of @ such that Xe S 'I'p n 


Proof. The set 4'p fl is compact, forward invariant, and convex, since the both sets 'Pp,'^ are as such. 
Hence, the statement of the lemma follows directly from the application of the Brouwer’s fixed point theorem 
31j on the associated flow. □ 

We are ready now to prove the second item of Theorem ITT] 


Proof of Theorem [TJf 2). By Lemma 1181 there exists an equilibrium in 5'p. Since it is assumed that 
there exists an isolated equilibrium in interior, Theorem 1101 implies that is not GS. □ 


Before concluding the proof, a simple lemma is stated and proved: 

Lemma 19. Let Xe be an equilibrium of ([3]). Let P be a set of species that correspond to {1, ..,n}\ supp(xe). 
Then, P is a siphon. 

Proof. Assume that P is not a siphon, then there exists some Xi G P and Rj € Sf such that Xi is a product 
of Rj and Rj 7 ^ A(P). At the given equilibrium, all negative terms in the expression of Xi vanish since 
Xei = 0. Since Xi is not a reactant in R_, this implies /3ij > 0,ay = 0 then Rj{x) has a strictly positive 
coefficient which implies Xi >0; a contradiction. □ 

Hence, we are ready to conclude the proof of Theorem 1111 

Proof of Theorem \ll\ -3). By Lemma [T31 there exists an equilibrium x* G p such that ri?(x*) = 0. 
Since dim(kerr) = 1, this implies that R{x*) = tv for some t > 0. Consider the case t = 0. This implies 
R{x*) = 0. Then, P G P := {1,.., n}\supp(x*). P is a siphon by Lemma [T9l and since P C P it is a critical 
deadlock. However, by Theorem ITTl-ll. gVa^b is not GS. If t > 0, this implies that P = 0; a contradiction. □ 


Proof of Theorem B3I Write the logarithmic norm expression using dini): 


hc{Ji{0) = limsup - 
ft—>-0+ ^ 


ax 



The expression above includes the induced matrix norm. However, it is not easy to find a explicit expression 
of the induced matrix norm with respect to non-square weighting. Therefore, the definition of the induced 
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matrix norm is used as follows: 




= sup 

C IICClloo-1 


= sup 

(♦) 

= sup 
IIC^IIoo-1 




Ci + 

£^1 
s 

C^ + hJ2peA^Ci 


e=i 


< sup 
liciiioo=i 


e=i 


I + hJ2peA^ 
I + hY,Pe^‘ 


relic 


i=l 


where (*) is by ® and (X) is by (EZD. Therefore, the expression of the logarithmic norm above can be 
written as: 

PciMO) < Poo i^peA^j <^PiPooiA^) = 0, (28) 

Vfci / e=i 

where the inequalities follow by the subadditivity of the logarithmic norm and m- 

Note since C has nonempty kernel space, then 11(7^1100 is a semi-norm. However, Theorem 1131 requires 
a norm. This can be remedied by studying the system in directions orthogonal to kerT by defining a 
transformation of coordinates using a matrix Ti = [Ti, ui,.., where {ui,.., v^-ur-} is a basis of kerT, 

Vr = rank{T), and Ti is chosen so that Ti is invertible. Defining ^ = T^, then the first Vr coordinates of ^ are 
decoupled from the rest. Hence, inequality (E51) can be established similarly for the reduced Jacobian which is 
the upper right Vr x Vr block of TiJi{^)T^^. The norm in the reduced subspace is ||.||(=, : z HCzHoo, ^ € M*", 
where C is ^ ^ ’^r defined as the nonzero columns of CT-^^. The equations = MC are equivalent to 

(Crj“^)(riej 7 )^rj“^) = A^{CT^^). Therefore, everything goes through in the reduced space, and the same 
upper bound is valid. 

The same argument can be replicated for to prove the second item in the theorem. This is accomplished 
by utilizing the alternative representation (|8]), the rank-one decomposition of the Jacobian ([9]), and noting 
that conditions (ED can be written as: BTicJ = A^B + for some matrices The reduced 

space argument can be carried out also by using a transformation matrix T 2 = [T 2 , . ■ 
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